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Abstract. Middle East Respiratory Syndrome Coronavirus (MERS-CorV), was discovered in humans with lower respiratory tract 
infection, causes a range of illnesses in humans, from the common cold to the Severe Acute Respiratory Syndrome (SARS). 
Scientists give much attention to study the CorV infection among groups and travelers. In this paper, we utilize a mathematical 
model governed by a system of differential equations, which incorporate target cell limitation and the innate interferon response, 
investigate the innate and adaptive immune responses to primary CorV infection in an individual. We also investigate the sensitivity 
analysis of the model to determine the most sensitive parameters and informative subintervals. This study may promote clearance 
of virus and host recovery from infection. 
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Introduction 


Infectious diseases remain a major cause of death, disability, and social and economic disorder for millions of people 
throughout the world. Prevention and treatment strategies for infectious diseases are derived thorough understanding 
of the complex interactions between specific viral or bacterial pathogens and the human (or animal) host. Recently, 
there are increasing attempts to use mathematical modelling tools in biology and medicine. A particular interest to 
mathematical models has been established in theoretical studies of the immune system and infectious diseases in the 
molecular level. Mathematicians and immunologists have begun to work together to create models that attempt to 
predict the progression of disease in an individual, using differential models; See [1, 2, 3, 4, 5, 6]. 

Herein, we give attention to human Coronavirus (CorV) infections, which first appeared in April 2012 in the 
Middle East. Since then, 536 laboratory-confirmed cases of human infection with Middle East respiratory syndrome 
coronavirus (MERS-CorV) have been reported to WHO, including 145 deaths [7, 8, 9]. To date, the affected countries 
in the Middle East include Jordan, Kuwait, Oman, Qatar, Saudi Arabia (KSA), United Arab Emirates (UAE) and 
Yemen; in Africa: Egypt and Tunisia; in Europe: France, Germany, Greece, Italy and the United Kingdom; in Asia: 
Malaysia and Philippines; and in North America: the United States of America (USA). All of the cases recently 
reported outside the Middle East (Egypt, Greece, Malaysia, the Philippines and the USA) recently travelled from 
countries inside of the Middle East (KSA or UAE). Overall, 65.6% of cases are male and the median age is 49 years 
old (range 9 months-94 years old). Since the last update of 27 March 2015, 330 laboratory-confirmed cases, including 
59 deaths, were reported to WHO. These cases have not reported contacts with other laboratory-confirmed cases, and 
some have reported contacts with animals, including camels. Although camels are suspected to be the primary source 
of infection for humans, the routes of direct or indirect transmission remain unknown and investigations are ongoing 
[10, 11]. 

Mathematical models have proven to be useful tools in the analysis of viral infections. For example, ordinary and 
partial differential equations have long played important roles in bioscience, and they will no doubt continue to serve 
as indispensable tools in future investigations. Although studying transmission of the diseases among groups (animals 
and humans) using the epidemic (or endemic/ pandemic) SIR or SIS models have a long and well-established history 
of research, deterministic mathematical modeling of infectious diseases on the molecular level is still a relatively new 
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field (see, e.g., [12, 13, 14, 15, 16, 17]). The concern in this paper is to provide a system of ODEs model to describe 
the dynamic process of the interactions between the immune system and CorV infections in an organism. Qualitative 
behaviour of the model is investigated by sensitivity analysis of the parameters to model states. 


The organization of this paper is as follows.. In Section 2, we introduce a general background about the interaction 
of immune system with CorV infection. In Section 3, we propose a mathematical model based on a system of ODEs. 
In Section 4, we briefly introduce sensitivity analysis of the model. 


General Background 


Coronavirus (CorV) attacks the host respiratory tract mucosa, interacts with healthy epithelial cells and infects them 
by binding to cell surface receptors via one of the major surface glycoproteins, HA. The virus replicates in infected 
cells and several hours after cellular infection, newly synthesized virus particles are released by the action of another 
major glycoprotein, NA. The response of the host to CorV infection involves a cascade of events mediated by several 
effector cells and molecules that neutralize free virus, kill infected cells and limit the spread of viral particles by 
increasing healthy cell resistance to infection [18, 19]. 


Antigen presenting cells (APC) are essential in the induction and amplication of the human immune response. 
Exogenous viral antigens, which comprise inactive viral particles, intact viruses and apoptotic, infected cells, are 
taken up by APC through endocytosis and provide a potential source of peptides that could bind to MHC class I or 
II molecules in the APC. The role of the APC is to stimulate both innate and adaptive immunity. As the rst line of 
defense, APC and infected cells stimulate the innate immunity by secreting interferon a and 6 (IFN) molecules which 
interact with healthy cells and convert them to an infection resistant state, thereby preventing the virus from spreading 
efciently and allowing the adaptive immune response enough time to develop and eliminate the virus.Another role of 
IFN is to stimulate symptoms such as fever which occurs in the early stages of infection. IFN levels rise rapidly after 
infection and correlate directly with the degree of viral replication in ferrets, mice and humans. Magnitude of the fever 
correlates strongly with the level of virus shedding in humans and animals [20]. 


As a second line of defense, APC stimulate the cellular component of innate immunity which consists of effector 
cells (cytotoxic T cells (CTL) or natural killer cells (NK)) that destroy infected cells before they can release a mature 
virus. Activated T cells produce various factors which are extremely important for the kinetics of the CorV infection: 
helper T cells secrete IL-2 and other lymphocytes and CTL produce IFN-y, which increases the expression of MHC 
antigens acting to enhance virus-infected cell destruction. The peptide-class I MHC complexes presented on the 
infected cells are recognized by class-I MHC-restricted CD8+ memory T-cells (Th1 cells), which destroy the infected 
cells. The specicity of memory T cells is directed against viral internal proteins; NP is the strongest of these antigens. 
Since structure of these antigens is conserved within the type of virus, Th1 cells against these antigens are cross- 
reactive within the type of virus [21]. 


Finally, APC stimulate adaptive immunity by activating the proliferation of virus-specic plasma cells which 
produce antibodies (Abs) that bind with CorV and render it ineffective. HA and NA are taken up in an endocytic 
vesicle pathway of the APC and are degraded; the peptides of these antigens are loaded on class- MHC molecules 
and then expressed on the APC (Tulp et al., 1994). The peptide class II MHC complexes are recognized by class-II 
MHC-restricted CD4+ T cells (Th2 cells). Th2 cell stimulation by antigen recognition results in the produc- tion 
of specic Abs to HA and NA molecules. Anti-HA Abs neutralize the infectivity of the virus, whereas anti-NA Abs 
prevent the release of viruses from infected cells. Thus, anti-HA Abs are primarily responsible for preventing infection, 
while anti-NA Abs and CTL specic for viral core proteins are responsible for reducing viral spread and thereby for 
accelerating the recovery from the virus. 


The respiratory tract mucosa is not only the site of infection by inuenza viruses but also the site of defense 
against viral infection in the host. The recovery process after primary infection involves two phases: an early phase 
(days 57), characterized by a rapid decrease in virus titer via killing of the virus- infected epithelial cells by MHC 
class I restricted CD8+ CTLs, which appear with a peak at day 7 is cellular response dependent, while a late phase 
(day 7 onwards), characterized by a more protracted decrease that ultimately results in clearance, depends on the 
adaptive response. Since u symptoms emerge within a few days of inoculation, acquired immunity appearing after 
the rst week of infection cannot prevent the onset of respiratory symptoms. Therefore, effective immunity must be 
induced in advance by natural infection or vaccination in order to prevent disease [8, 22, 23]. 
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FIGURE 1. Schematic representation of interactions included in model (1). The CorV is shown as red hexagon, the four different 
cell types are shown in cyan. Components of adaptive immunity are shown in orange, cellular component of innate immunity 
in purple, and interferon component in green. Up-regulation is represented by lines terminated with arrows and inhibition by 
lines terminated with bars. The inter-conversion of cell types is indicated by dashed arrows. Dashed ovals represent details of the 
mechanism ignored in this model. Homeostatic maintenance of effector and plasma cell populations are indicated by self-regulating 
loops. 


The Model 


In this paper, we construct a mathematical model of the dynamics of CorV infection and the human immune response 
to such infection. We focus on three important components of the immune response: the interferon and cellular com- 
ponents of innate immunity and the adaptive immunity, all of which have the same goal of limiting the concentration 
of the virus and the damage to the system, but which achieve this goal using different strategies: interferon immunity 
by removing the sub- strate that virus needs for reproduction (1.e., the healthy cells), cellular immunity by removing 
the source of new viruses (1.e., the infected cells), and adaptive immunity by lowering the effective concentration of 
the virus. 

The model of human immune response against CorV infection we consider is a simplied model of population- 
dynamics type which consists of the following interactions (see Fig. 1). The epithelial cells of the respiratory tract 
are assumed to be in one of four possible states: healthy (H), infected (1), dead (D), or resistant (R) to infection. The 
total number of epithelial cells G.e., H+I+D+R) is assumed constant. The virus particles (V) interact with healthy 
cells and infect them. Infected cells release new virus particles upon their death. Proliferation of healthy cells causes 
regeneration and decrease in the proportion of dead cells. Dead cells stimulate the activation of APC (M). 

The interactions are based on clonal selection theory, massaction kinetics, characteristics of interactions and the 
birthdeath balances of populations of cells and molecules. Given the variables of table 1, we assume the following 
model which is based on th emodel is similar on the model that was proposed by Bocharov and Romanyukha [24]. 


VO = yl) -yvvaS OAOVO -yv HOVE) -ay V(t) — we (1) 
A(t) = bypD (A) + RW) + arR@® — yuv VOHA) — bur FOH(), (2) 

rO = yav VOHA) -b EOI — al), (3) 
M(t) = (bypD(t) + bmv VAA - M()) - amM(t), (4) 
Ft) = bpM(t)+ crI(t)-cryHOF(t) -ar Fo), (5) 
R(t) = byrF(t)H(t) — arR(t), (6) 
Et) = bemMQE() — be TOE) — ag — EO), (7) 
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P(t) = bpmM(t)P(t) -apC - PO), (8) 
At) = bgA(t) — YavS AVOA) — a4 A(t), (9) 
S'HA = rP -S(d). (10) 
Here D(t) is the proportion of dead cells which is given by 
D(t) =1- H(t) - R(t) - 1M. (11) 
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FIGURE 2. Shows, for particular values of parameters, the numerical simulation of the model (1-11). 


First Eq. of (1) of the system describes the rate of change of virus concentration V. It expresses the production 
rate of a viral particle by infected cells, rate of neutralization of CorV by specific antibodies, the rate of adsorption 
of viral particles by healthy cells, and the natural decay of viral particles. The viral particles are also removed from 
the respiratory tract by nonspecific mechanisms. The nonspecific mucociliary removal of virions supported by cough 
and other mechanisms is described by the term ayı V/(1 + ay2V), which saturates with increasing V as the available 
capacity of these mechanisms is exhausted. Note that the lethal damage of an infected cell by the effector cells does 
not cause any release of infective CorV and hence in Eq. (1) there is no term of the form El. 

Second Eq. of (2) determines the time rate of change of healthy cells H. During recovery, new healthy cells are 
generated as a result of proliferation of both healthy and resistant cells (the offspring of resistant cells lose resistance) 
and hence the proliferation term is proportional to (H + R), and to D (in a logistic fashion) since regeneration can only 
occur in the presence of damage. Resistant cells R gradually lose their resistance to infection and return into their 
initial sensitive state (healthy state) (Joklik, 1985), which is characterized by the term arR. The term gyyVA is the 
loss of healthy cells due to infection and the term byr FH characterizes transition of the healthy cells into resistant 
state. Third Eq. of (3) characterizes the time rate of change of infected cells 7. The infection of healthy cells by virions 
is described in the term yyy VH. The term a;/ indicates the natural death of infected cells during which new virus 
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FIGURE 3. Displays cumulative proportions of types of respiratory epithelial cells. For particular values of parameters It shows, 
at any given time, the proportions of dead cells, infected cells, resistant cells, and healthy cells. 


TABLE 1. Variables of Model (1) [3]. 


| Variable | Description 


VÆ 
A(t) 
I(t) 

M(t) 
F(t) 
R(t) 
E(t) 
P(t) 
A(t) 
S(t) 


Viral load per epithelial cell 

Proportion of healthy cells 

Proportion of infected cells 

Activated antigen presenting cells per homeostatic level 
Interferon per homeostatic level of macrophages 
Proportion of resistant cells 

Effector cells per homeostatic level 

Plasma cells per homeostatic level 

Antibodies per homeostatic level 

Antigenic distance 








Sensitivity function, d Vð W 








Sensitivity function, ð V/@ Wu 
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FIGURE 4. shows the sensitivity V(t) due to small perturbations on the parameters yy, and yyy. The state variable V(t) is insensi- 
tive the beginning then sensitive to small changes in the parameters. 


particles are produced. The term bzg EI characterizes the destruction of infected cells by effector cells (CTL and NK) 
during which no new virus is produced. 

Fourth Eq. of (4) establishes that the time rate of increase of activated APC (M) is proportional to the amount 
of the virus and the amount of dead cells. The natural decay of activated state of APC is represented by the last term 
in that equation. The sixth Eq. describes the time rate of change of interferon œ and 6 (F) which depends on the 
production rate of F by APC and by infected cells, on the rate of F binding healthy cells, as well as on the nonspecific 
decay of F. The fifth Eq. shows that resistant cells R are induced from healthy cells (byr FH) and convert back to 
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healthy cells (arR) with finite lifetime. 

Seventh Eq. of (7) characterizes the rate of change of effector cells E concentration and takes into account the 
production rate of effector cells stimulated by APC (first term) and the destruction rate of infected cells by effector 
cells (second term, bg; IE). The terms ag(1 — E) and ap(1 — P) and (8) are approximate expressions for homeostatic 
maintenance of the levels of active effectors and plasma cells, reflecting the observation that the healthy body tends to 
maintain their concentrations within narrow bounds. In a healthy state the effectors and plasma cells are naturally lo- 
cated in lymph nodes and blood, and migrate into the infected tissue upon activation. Both the activation and migration 
of those cells to the infected tissue are assumed to be much faster than their proliferation and hence are not explicitly 
accounted for. The first term in eight Eq. of (8) characterizes the activation process of plasma cells stimulated by APC. 

Ninth Eq. of (9) stands for the time rate of change of the concentration of antibodies A describing the production 
rate of A by plasma cells (first term), the neutralization rate of free viral particles by specific antibodies (second 
term) and the natural decay rate of A (last term). The variable S in our model represents the compatibility between 
antibodies and the virus strain in an individual and ranges from O (no compatibility) to 1 (maximal compatibility) 
and can be interpreted as a measure of binding affinity of the antibody and the virus. The immune memory of the 
host is described by the initial value S (0) of S. During the course of the disease, S increases as plasma cells produce 
antibodies increasingly compatible with viral antigens. The rate of increase of S is approximated by the term rP(1 — S) 
which accounts for two natural observations: (1) the increase in S is stimulated by plasma cells and (11) S cannot 
increase beyond 1. By adjusting the time evolution of S we may observe how the course of the disease depends on the 
evolution of antigenic distance. 

We note that a number of assumptions in the model are strong simplifications of our knowledge of immune 
physiology. The populations of cells and virus are assumed to be uniformly distributed over the epithelial layer at 
all times. It is also assumed that time rate of change of any model variable is determined by the present value of all 
variables. Some of the variables do not have uniquely identifiable biological counterparts. For example, there is no 
single biological entity or marker that represents the APC in our model, which are assumed to provide both antigen 
presenting and JFN producing functions. We have also omitted intermediate steps in the pathways: for example, we 
do not account for the intermediate steps in the production of effector cells and plasma cells such as Thl and Th2 
helper cells and B-cells. We do not consider time delays in the reproduction of cellular components. 

Figure 2 shows the numerical simulation of the model, or particular values of parameters. While Figure 3 displays 
cumulative proportions of types of respiratory epithelial cells: at any given time, the proportions of dead cells, infected 
cells, resistant cells, and healthy cells. 


Sensitivity Analysis 


The goals of sensitivity analysis with respect to random perturbations of the model parameters are (i) to show how 
robust of the the infection model is in relation to perturbed parameter values, (ii) to explore to which parameters the 
system is more sensitive to understand key processes and immune system mechanisms. Figure 4 shows the sensitivity 
V(t) due to small perturbations on the parameters yy, and yyy. The state variable V(t) is insensitive the beginning then 
sensitive to small changes in the parameters; See [22, 25, 26]. 


Conclusion 


We presented a mathematical model, based on a system of ODEs, of the human immune response to the CorV infection 
in individual hosts which includes innate and adaptive immunity, and analyzed its behavior. Such a model could be 
used to explore in more detail individual determinants of symptoms and behavior of clinical relevance, especially in 
large-scale simulations of disease spread and containment. Simulation and sensitivity analysis of this model suggest 
that for majority of possible parameter values and initial conditions the course of the disease falls into one of three 
categories: asymptomatic disease, typical disease, and severe disease. The model can be extended to consider time 
delays in the reproduction of cellular components. 
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